{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Dataflowr](https://raw.githubusercontent.com/dataflowr/website/master/_assets/dataflowr_logo.png)](https://dataflowr.github.io/website/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "GPfXTpAXB6MI"
   },
   "source": [
    "# Predicting Engine Failure with RNN\n",
    "\n",
    "In this practicals, the goal is to predict the failure of an engine. The training dataset is made of time series obtained from several sensors on the engine until failure. The test dataset is made of the start of these time series and a failure date.\n",
    "\n",
    "We will build a simple RNN taking as input the multi-dimensional time serie characterizing the engine and learn its parameters to predict the time of failure at each instant. At the start, the best prediction without any input data should be the average of the failure times in the dataset and as more and more data is fed in the RNN, it should give a better and better estimate.\n",
    "\n",
    "The dataset is provided by [NASA](https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/#turbofan), see also [Kaggle](https://www.kaggle.com/datasets/suriyachayatummagoon/cmapssdata?select=Damage+Propagation+Modeling.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "egcRZderB6MK"
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.special import gamma\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ZZD9lpogB6MP"
   },
   "outputs": [],
   "source": [
    "torch.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nCZo-JIUB6MT"
   },
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "print('Using gpu: %s ' % torch.cuda.is_available())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_tMC18fSB6MW"
   },
   "source": [
    "## 1. Downloading the data\n",
    "\n",
    "This need to be done only once!\n",
    "\n",
    "You can find the data on the website of the [NASA](https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/#turbofan) or on [Kaggle](https://www.kaggle.com/datasets/suriyachayatummagoon/cmapssdata?select=Damage+Propagation+Modeling.pdf) or on my website:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6IUyhYERB6MX"
   },
   "outputs": [],
   "source": [
    "%mkdir data\n",
    "%cd data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "xX0tNgXvB6Ma"
   },
   "outputs": [],
   "source": [
    "!wget 'https://www.di.ens.fr/~lelarge/CMAPSSData.zip'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4-R0S-6AB6Mc"
   },
   "outputs": [],
   "source": [
    "!unzip CMAPSSData.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "bDsmGAAKB6Mf"
   },
   "outputs": [],
   "source": [
    "%cd .."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jw7cICW_B6Mi"
   },
   "source": [
    "## 2. Loading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "LgNZ-5pzB6Mj"
   },
   "outputs": [],
   "source": [
    "def get_CMAPSSData(nb_file):\n",
    "    # get data from file and pre process it (normalization and convert to pandas)\n",
    "    dataset_train = pd.read_csv('./data/train_FD00{}.txt'.format(nb_file),\n",
    "                                sep=' ', header=None).drop([26, 27], axis=1)\n",
    "    dataset_test = pd.read_csv('./data/test_FD00{}.txt'.format(nb_file),\n",
    "                               sep=' ', header=None).drop([26, 27], axis=1)\n",
    "    test_truth = pd.read_csv('./data/RUL_FD00{}.txt'.format(nb_file),\n",
    "                             sep=' ', header=None).drop([1], axis=1)\n",
    "    col_names = ['id', 'cycle', 'setting1', 'setting2', 'setting3', 's1', 's2', 's3', 's4', 's5', 's6', 's7', 's8',\n",
    "                 's9',\n",
    "                 's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's18', 's19', 's20', 's21']\n",
    "    dataset_train.columns = col_names\n",
    "    dataset_test.columns = col_names\n",
    "    test_truth.columns = ['more']\n",
    "    test_truth['id'] = test_truth.index + 1\n",
    "    rul = pd.DataFrame(dataset_test.groupby('id')['cycle'].max()).reset_index()\n",
    "    rul.columns = ['id', 'max']\n",
    "    test_truth['rtf'] = test_truth['more'] + rul['max']\n",
    "    test_truth.drop('more', axis=1, inplace=True)\n",
    "    dataset_test = dataset_test.merge(test_truth, on=['id'], how='left')\n",
    "    dataset_test['ttf'] = dataset_test['rtf'] - dataset_test['cycle']\n",
    "    dataset_test.drop('rtf', axis=1, inplace=True)\n",
    "    dataset_train['ttf'] = dataset_train.groupby(['id'])['cycle'].transform(max) - dataset_train['cycle']\n",
    "    features_col_name = ['setting1', 'setting2', 'setting3', 's1', 's2', 's3', 's4', 's5', 's6', 's7', 's8',\n",
    "                         's9', 's10', 's11',\n",
    "                         's12', 's13', 's14', 's15', 's16', 's17', 's18', 's19', 's20', 's21']\n",
    "    target_col_name = 'ttf'\n",
    "    relevant_features_col_name = []\n",
    "    for col in features_col_name:\n",
    "        if not (len(dataset_train[col].unique()) == 1):\n",
    "            relevant_features_col_name.append(col)\n",
    "    sc = MinMaxScaler()\n",
    "    dataset_train[features_col_name] = sc.fit_transform(dataset_train[features_col_name])\n",
    "    dataset_test[features_col_name] = sc.transform(dataset_test[features_col_name])\n",
    "    return dataset_train, dataset_test, relevant_features_col_name, target_col_name\n",
    "\n",
    "\n",
    "def to_lists_of_tensors(dataset, features_col_name, target_col_name):\n",
    "    # take pandas df and convert it to list of tensors (for pytorch)\n",
    "    X, y = [], []\n",
    "    nb_sequences = max(dataset['id'])\n",
    "    for i in range(1, nb_sequences + 1):\n",
    "        df_zeros = dataset.loc[dataset['id'] == i]\n",
    "        df_one_x = df_zeros[features_col_name]\n",
    "        df_one_y = df_zeros[target_col_name]\n",
    "        X.append(torch.from_numpy(np.expand_dims(df_one_x.values, 1)).type(torch.FloatTensor))\n",
    "        y.append(torch.from_numpy(df_one_y.values).type(torch.FloatTensor))\n",
    "    return X, y\n",
    "\n",
    "\n",
    "def convert_train_and_test_to_appropriate_format(dataset_train, dataset_test, features_col_name, target_col_name):\n",
    "    # take 2 datasets (train and test and covert them to lists of tensors)\n",
    "    X_train, y_train = to_lists_of_tensors(dataset_train, features_col_name, target_col_name)\n",
    "    X_test, y_test = to_lists_of_tensors(dataset_test, features_col_name, target_col_name)\n",
    "    return X_train, y_train, X_test, y_test\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "4IJ-VFPXB6Ml"
   },
   "outputs": [],
   "source": [
    "%pycat ./data/readme.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "X-xRkXcQB6Mn"
   },
   "outputs": [],
   "source": [
    "dataset_train, dataset_test, features_col_name, target_col_name = get_CMAPSSData(1)\n",
    "X_train, y_train, X_test, y_test = convert_train_and_test_to_appropriate_format(dataset_train, dataset_test,\n",
    "                                                                                    features_col_name, target_col_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "uJoXsHBwB6Mq"
   },
   "outputs": [],
   "source": [
    "dataset_train.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Bs3X-PShB6Mt"
   },
   "source": [
    "Here I have done the minimal preprocessing of the data using the [`MinMaxScaler`](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.MinMaxScaler.html) of sklearn to scale each feature in (0,1).\n",
    "\n",
    "`X_train` is a list where each element is of shape (length_of_sequence,1,number_of_sensors) where the second dimension with value 1 corresponds to the batch size. As in the course, we will not proceed sequences by batches but one after the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "02WrnXXEB6Mt"
   },
   "outputs": [],
   "source": [
    "X_train[0].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ayDGRokNB6Mw"
   },
   "source": [
    "## 3. WTTE-RNN model\n",
    "\n",
    "Here, we follow an approach inspired from this [blog](https://ragulpr.github.io/2016/12/22/WTTE-RNN-Hackless-churn-modeling/#wtte-rnn-produces-risk-embeddings).\n",
    "\n",
    "You first need to define a GRU (or LSTM) that will take as input a sequence of shape (length_of_sequence,1,number_of_sensors) and output a sequence of shape (length_of_sequence,2) obtained by passing the output of the GRU through a linear layer. As we want positive number, you will take the exponent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "61SxB3gYB6Mw"
   },
   "outputs": [],
   "source": [
    "class GRUnet(nn.Module):\n",
    "    def __init__(self, dim_input, num_layers, dim_hidden, dim_output=2):\n",
    "        super(GRUnet, self).__init__()\n",
    "        #\n",
    "        # your code here\n",
    "        #\n",
    "        \n",
    "\n",
    "    def forward(self, x):\n",
    "        #\n",
    "        # your code here\n",
    "        #\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ezsnVZjLB6Mz"
   },
   "source": [
    "Test your network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NQAWWJ_kB6Mz"
   },
   "outputs": [],
   "source": [
    "model = GRUnet(dim_input=len(features_col_name), num_layers=3, dim_hidden=50)\n",
    "model = model.to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DLlNMl69B6M2"
   },
   "outputs": [],
   "source": [
    "output = model(X_train[0].to(device))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nJX0-hnjB6M4"
   },
   "outputs": [],
   "source": [
    "output.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AhHb7vleB6M7"
   },
   "source": [
    "In order to learn the parameters of your RNN, you need to specify a loss and we will here follow a standard approach in reliability theory: we model the failure time as a [Weibull random variable](http://reliawiki.org/index.php/The_Weibull_Distribution)\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(X>t) = \\exp\\left( \\frac{t}{\\eta}\\right)^{\\beta},\n",
    "$$\n",
    "where $\\eta$ is the scale parameter and $\\beta$ is the shape parameter.\n",
    "\n",
    "Note that we have for the mean of a Weibull distribution:\n",
    "$$\n",
    "\\mathbb{E}[X] = \\eta \\Gamma(1+1/\\beta),\n",
    "$$\n",
    "where $\\Gamma$ is the [Gamma function](https://en.wikipedia.org/wiki/Gamma_function).\n",
    "\n",
    "In our case, we will interpret the 2 outputs of the RNN as estimates for the parameters $\\eta$ and $\\beta$. In order to design a loss, we compute the log-likelihood:\n",
    "\\begin{eqnarray*}\n",
    "\\log f(t) &=& \\log\\left( \\frac{\\beta}{\\eta}\\right) +(\\beta -1)\\log\\left(\\frac{t}{\\eta}\\right) -\\left(\\frac{t}\n",
    "{\\eta} \\right)^{\\beta}\\\\\n",
    "&=& \\log \\beta +\\beta \\log\\left(\\frac{t}{\\eta}\\right) -\\log t-\\left(\\frac{t}\n",
    "{\\eta} \\right)^{\\beta}\n",
    "\\end{eqnarray*}\n",
    "\n",
    "Define a loss function corresponding to the negative log-likelihood (add a small parameter $\\epsilon$ to $t$ in order not to compute $\\log 0$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "mm7PPLLZB6M8"
   },
   "outputs": [],
   "source": [
    "class weibull_loss(nn.Module):\n",
    "    def __init__(self):\n",
    "        super(weibull_loss, self).__init__()\n",
    "        self.epsilon = 1e-6\n",
    "\n",
    "    def forward(self, output, y):\n",
    "        #\n",
    "        # your code here\n",
    "        #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "AJEqkoBsB6M-"
   },
   "source": [
    "Test your loss function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "QkyZMThqB6M-"
   },
   "outputs": [],
   "source": [
    "loss_fn = weibull_loss()\n",
    "loss_fn(output.squeeze(),y_train[0].to(device))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VxdK_B1vB6NB"
   },
   "source": [
    "## 4. Training your model\n",
    "\n",
    "Code your taining and testing loops.\n",
    "\n",
    "You might want to use a scheduler like `torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, factor=0.5, patience=1, verbose='True',threshold=0.001)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "nR1GA2KzB6NC"
   },
   "outputs": [],
   "source": [
    "def train_epoch(X_train, y_train, model, loss_fn, optimizer, device):\n",
    "    # train the model through the whole training dataset for one epoch\n",
    "    # return the corresponding loss on the epoch\n",
    "    #\n",
    "    # your code here\n",
    "    #\n",
    "\n",
    "\n",
    "def test_epoch(X_test, y_test, model, loss_fn, device):\n",
    "    # evaluate the model through the whole testing dataset\n",
    "    # return the corresponding loss\n",
    "    #\n",
    "    # your code here\n",
    "    #\n",
    "\n",
    "\n",
    "def fit(model, X_train, y_train, X_test, y_test, optimizer, loss_fn, nb_epochs, device):\n",
    "    # fit the model by training it nb_epochs times\n",
    "    # you might want to use a scheduler\n",
    "    #\n",
    "    # your code here\n",
    "    #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "lBixgIwBB6NE"
   },
   "outputs": [],
   "source": [
    "model = GRUnet(dim_input=len(features_col_name), num_layers=3, dim_hidden=50,dim_output=2)\n",
    "model = model.to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "lwYUv3t5B6NG"
   },
   "outputs": [],
   "source": [
    "optimizer = torch.optim.Adam(model.parameters(), lr=0.01)\n",
    "loss_fn = weibull_loss()\n",
    "nb_epochs = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1J8VDkBTB6NI"
   },
   "outputs": [],
   "source": [
    "model, train_loss_t, test_loss_t = fit(model, X_train, y_train, X_test, y_test, optimizer, loss_fn, nb_epochs,device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "HL--b9nVB6NK"
   },
   "outputs": [],
   "source": [
    "def plot_losses(train_loss_t, test_loss_t):\n",
    "    nb_epochs = len(train_loss_t)\n",
    "    plt.plot(range(nb_epochs), train_loss_t, color='orange', label='Loss on the training set')\n",
    "    plt.plot(range(nb_epochs), test_loss_t, color='green', label='Loss on the testing set')\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "op0o9tx4B6NN"
   },
   "outputs": [],
   "source": [
    "plot_losses(train_loss_t, test_loss_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eFZ9dY4HB6NQ"
   },
   "source": [
    "## 5. Looking at your results\n",
    "\n",
    "To compute a baseline, I am computing the average of all failure times in the train dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FVtYbGjkB6NQ"
   },
   "outputs": [],
   "source": [
    "max_val = np.zeros(len(y_train))\n",
    "for i,y in enumerate(y_train):\n",
    "    max_val[i] = y[0].item()\n",
    "baseline = np.mean(max_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xYqaD1loB6NS"
   },
   "source": [
    "Here I am computing all the predictions made by my model on the test set and exporting these into `numpy`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "F_UttARgB6NS"
   },
   "outputs": [],
   "source": [
    "def compute_np(model,X_test, y_test,baseline=baseline,device=device,max_size=303):\n",
    "    n_test = len(X_test)\n",
    "    all_pred = np.empty((n_test,max_size,2))\n",
    "    all_y = np.empty((n_test,max_size))\n",
    "    base_pred = np.empty((n_test,max_size))\n",
    "    all_pred[:] = np.NaN\n",
    "    all_y[:] = np.NaN\n",
    "    base_pred[:] = np.NaN\n",
    "    list_npred = []\n",
    "    for k in range(n_test):\n",
    "        pred = model(X_test[k].to(device))\n",
    "        pred_np = pred.cpu().detach().numpy()\n",
    "        n_pred = pred_np.shape[0]\n",
    "        list_npred.append(n_pred)\n",
    "        all_pred[k,:n_pred,:] = pred_np\n",
    "        all_y[k,:n_pred] = y_test[k].numpy()\n",
    "        base_pred[k,:n_pred] = baseline - range(n_pred)\n",
    "    return all_pred, all_y, base_pred, list_npred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6XWfOw9IB6NU"
   },
   "outputs": [],
   "source": [
    "all_pred, all_y, base_pred, list_npred = compute_np(model, X_test,y_test)\n",
    "pred_fail = all_pred[:,:,0]*gamma(1+1/all_pred[:,:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "3IdjZMAZB6NW"
   },
   "source": [
    "On the test set, we have only access to the start of the sequence and need to predict the failure time. For a given engine, you can compare the predictions made by the model, the baseline and the true value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CdkmllHNB6NX"
   },
   "outputs": [],
   "source": [
    "k = 50\n",
    "plt.plot(pred_fail[k,:],label='predicted')\n",
    "plt.plot(all_y[k],label='true')\n",
    "plt.plot(base_pred[k],label='baseline')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oeFL3HnQB6NZ"
   },
   "source": [
    "To get a measure of perfomance, we compute the RMSE: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "NKMLIflYB6Na"
   },
   "outputs": [],
   "source": [
    "def RMSE(pred_fail, all_y):\n",
    "    return np.sqrt((pred_fail-all_y)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "yAnsK8GxB6Nc"
   },
   "outputs": [],
   "source": [
    "res= RMSE(pred_fail,all_y)\n",
    "res_base = RMSE(base_pred,all_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mrToJiZTB6Ne"
   },
   "source": [
    "The RMSE error is the distance between the esimation and the true line above. It is constant for the baseline and should decrease as we get more and more data with our model. Here is an example on the particular exmaple above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "VaCpHTDvB6Ne"
   },
   "outputs": [],
   "source": [
    "plt.plot(res[k])\n",
    "plt.plot(res_base[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "9sta7kgwB6Ng"
   },
   "source": [
    "Below, I am averaging the RMSE over the dataset keeping the time axis (note that each point is an average but not with the same number of samples).\n",
    "\n",
    "We see that the RMSE of the baseline is very bad for long sequences. This should be expected as these long sequences corresponds to healthy engines!\n",
    "\n",
    "To the contrary, our model get a decreasing RMSE as a function of the length of the input sequence.\n",
    "\n",
    "![](https://raw.githubusercontent.com/mlelarge/dataflowr/master/PlutonAI/img/rmse.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "7edZe3XAB6Nh"
   },
   "outputs": [],
   "source": [
    "plt.plot(np.nanmean(res,0), label = 'RMSE')\n",
    "plt.plot(np.nanmean(res_base,0), label ='RMSE baseline')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3MypycQKB6Nj"
   },
   "outputs": [],
   "source": [
    "np.nanmean(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "K7w8lxubB6Nk"
   },
   "outputs": [],
   "source": [
    "np.nanmean(res_base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "31ozUyp1B6Nm"
   },
   "outputs": [],
   "source": [
    "last_indices = list((~np.isnan(res)).sum(axis = 1) - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XvMJJAvtB6No"
   },
   "outputs": [],
   "source": [
    "np.mean([res[i,j] for i,j in enumerate(last_indices)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "MAF4CPtQB6Nq"
   },
   "outputs": [],
   "source": [
    "np.mean([res_base[i,j] for i,j in enumerate(last_indices)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JoFkoNORB6Ns"
   },
   "source": [
    "Above we see that we divided the RMSE by more than 2 compare to the baseline.\n",
    "\n",
    "Here we do a scatter plot of the predictions vs true values for the baseline method (closer to the diagonal in blue is better):\n",
    "\n",
    "![](https://raw.githubusercontent.com/mlelarge/dataflowr/master/PlutonAI/img/base_scatter.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "65JosTRLB6Ns"
   },
   "outputs": [],
   "source": [
    "plot = sns.jointplot(x=[all_y[i,j] for i,j in enumerate(last_indices)],y=[base_pred[i,j] for i,j in enumerate(last_indices)],dropna=True,kind=\"kde\", n_levels=30, color=\"g\");\n",
    "plot.ax_joint.plot([0,150], [0,150], 'b-', linewidth = 2);\n",
    "plot.set_axis_labels('true', 'predicted');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "JZ2KXoo8B6Nu"
   },
   "source": [
    "Now we do the same scatter plot with our model (colser to the diagonal in blue is better). We see a great improvement.\n",
    "\n",
    "![](https://raw.githubusercontent.com/mlelarge/dataflowr/master/PlutonAI/img/model_scatter.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "I1-qdjrkB6Nu"
   },
   "outputs": [],
   "source": [
    "plot = sns.jointplot(x=[all_y[i,j] for i,j in enumerate(last_indices)],y=[pred_fail[i,j] for i,j in enumerate(last_indices)],dropna=True,kind=\"kde\", n_levels=30, color=\"g\");\n",
    "plot.ax_joint.plot([0,150], [0,150], 'b-', linewidth = 2);\n",
    "plot.set_axis_labels('true', 'predicted');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "YuzDlaPiB6Nw"
   },
   "outputs": [],
   "source": [
    "plot = sns.jointplot(x=all_y,y=pred_fail,dropna=True,kind=\"kde\", space=0, color=\"g\");\n",
    "plot.ax_joint.plot([0,200], [0,200], 'b-', linewidth = 2);\n",
    "plot.set_axis_labels('true', 'predicted');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "C8AUCEBYB6Nx"
   },
   "source": [
    "[![Dataflowr](https://raw.githubusercontent.com/dataflowr/website/master/_assets/dataflowr_logo.png)](https://dataflowr.github.io/website/)"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "include_colab_link": false,
   "name": "11_Predictions_RNN_empty.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
